home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / liboctave / dColVector.cc < prev    next >
C/C++ Source or Header  |  1997-03-07  |  7KB  |  370 lines
  1. // ColumnVector manipulations.
  2. /*
  3.  
  4. Copyright (C) 1996 John W. Eaton
  5.  
  6. This file is part of Octave.
  7.  
  8. Octave is free software; you can redistribute it and/or modify it
  9. under the terms of the GNU General Public License as published by the
  10. Free Software Foundation; either version 2, or (at your option) any
  11. later version.
  12.  
  13. Octave is distributed in the hope that it will be useful, but WITHOUT
  14. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16. for more details.
  17.  
  18. You should have received a copy of the GNU General Public License
  19. along with Octave; see the file COPYING.  If not, write to the Free
  20. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  21.  
  22. */
  23.  
  24. #if defined (__GNUG__)
  25. #pragma implementation
  26. #endif
  27.  
  28. #ifdef HAVE_CONFIG_H
  29. #include <config.h>
  30. #endif
  31.  
  32. #include <iostream.h>
  33.  
  34. #include "f77-fcn.h"
  35. #include "lo-error.h"
  36. #include "mx-base.h"
  37. #include "mx-inlines.cc"
  38. #include "oct-cmplx.h"
  39.  
  40. // Fortran functions we call.
  41.  
  42. extern "C"
  43. {
  44.   int F77_FCN (dgemv, DGEMV) (const char*, const int&, const int&,
  45.                   const double&, const double*,
  46.                   const int&, const double*, const int&,
  47.                   const double&, double*, const int&,
  48.                   long);
  49. }
  50.  
  51. // Column Vector class.
  52.  
  53. bool
  54. ColumnVector::operator == (const ColumnVector& a) const
  55. {
  56.   int len = length ();
  57.   if (len != a.length ())
  58.     return 0;
  59.   return equal (data (), a.data (), len);
  60. }
  61.  
  62. bool
  63. ColumnVector::operator != (const ColumnVector& a) const
  64. {
  65.   return !(*this == a);
  66. }
  67.  
  68. ColumnVector&
  69. ColumnVector::insert (const ColumnVector& a, int r)
  70. {
  71.   int a_len = a.length ();
  72.   if (r < 0 || r + a_len > length ())
  73.     {
  74.       (*current_liboctave_error_handler) ("range error for insert");
  75.       return *this;
  76.     }
  77.  
  78.   for (int i = 0; i < a_len; i++)
  79.     elem (r+i) = a.elem (i);
  80.  
  81.   return *this;
  82. }
  83.  
  84. ColumnVector&
  85. ColumnVector::fill (double val)
  86. {
  87.   int len = length ();
  88.   if (len > 0)
  89.     for (int i = 0; i < len; i++)
  90.       elem (i) = val;
  91.   return *this;
  92. }
  93.  
  94. ColumnVector&
  95. ColumnVector::fill (double val, int r1, int r2)
  96. {
  97.   int len = length ();
  98.   if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
  99.     {
  100.       (*current_liboctave_error_handler) ("range error for fill");
  101.       return *this;
  102.     }
  103.  
  104.   if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  105.  
  106.   for (int i = r1; i <= r2; i++)
  107.     elem (i) = val;
  108.  
  109.   return *this;
  110. }
  111.  
  112. ColumnVector
  113. ColumnVector::stack (const ColumnVector& a) const
  114. {
  115.   int len = length ();
  116.   int nr_insert = len;
  117.   ColumnVector retval (len + a.length ());
  118.   retval.insert (*this, 0);
  119.   retval.insert (a, nr_insert);
  120.   return retval;
  121. }
  122.  
  123. RowVector
  124. ColumnVector::transpose (void) const
  125. {
  126.   return RowVector (*this);
  127. }
  128.  
  129. ColumnVector
  130. real (const ComplexColumnVector& a)
  131. {
  132.   int a_len = a.length ();
  133.   ColumnVector retval;
  134.   if (a_len > 0)
  135.     retval = ColumnVector (real_dup (a.data (), a_len), a_len);
  136.   return retval;
  137. }
  138.  
  139. ColumnVector
  140. imag (const ComplexColumnVector& a)
  141. {
  142.   int a_len = a.length ();
  143.   ColumnVector retval;
  144.   if (a_len > 0)
  145.     retval = ColumnVector (imag_dup (a.data (), a_len), a_len);
  146.   return retval;
  147. }
  148.  
  149. // resize is the destructive equivalent for this one
  150.  
  151. ColumnVector
  152. ColumnVector::extract (int r1, int r2) const
  153. {
  154.   if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  155.  
  156.   int new_r = r2 - r1 + 1;
  157.  
  158.   ColumnVector result (new_r);
  159.  
  160.   for (int i = 0; i < new_r; i++)
  161.     result.elem (i) = elem (r1+i);
  162.  
  163.   return result;
  164. }
  165.  
  166. // column vector by column vector -> column vector operations
  167.  
  168. ColumnVector&
  169. ColumnVector::operator += (const ColumnVector& a)
  170. {
  171.   int len = length ();
  172.  
  173.   int a_len = a.length ();
  174.  
  175.   if (len != a_len)
  176.     {
  177.       gripe_nonconformant ("operator +=", len, a_len);
  178.       return *this;
  179.     }
  180.  
  181.   if (len == 0)
  182.     return *this;
  183.  
  184.   double *d = fortran_vec (); // Ensures only one reference to my privates!
  185.  
  186.   add2 (d, a.data (), len);
  187.   return *this;
  188. }
  189.  
  190. ColumnVector&
  191. ColumnVector::operator -= (const ColumnVector& a)
  192. {
  193.   int len = length ();
  194.  
  195.   int a_len = a.length ();
  196.  
  197.   if (len != a_len)
  198.     {
  199.       gripe_nonconformant ("operator -=", len, a_len);
  200.       return *this;
  201.     }
  202.  
  203.   if (len == 0)
  204.     return *this;
  205.  
  206.   double *d = fortran_vec (); // Ensures only one reference to my privates!
  207.  
  208.   subtract2 (d, a.data (), len);
  209.   return *this;
  210. }
  211.  
  212. // matrix by column vector -> column vector operations
  213.  
  214. ColumnVector
  215. operator * (const Matrix& m, const ColumnVector& a)
  216. {
  217.   ColumnVector retval;
  218.  
  219.   int nr = m.rows ();
  220.   int nc = m.cols ();
  221.  
  222.   int a_len = a.length ();
  223.  
  224.   if (nc != a_len)
  225.     gripe_nonconformant ("operator *", nr, nc, a_len, 1);
  226.   else
  227.     {
  228.       if (nr == 0 || nc == 0)
  229.     retval.resize (nr, 0.0);
  230.       else
  231.     {
  232.       int ld = nr;
  233.  
  234.       retval.resize (nr);
  235.       double *y = retval.fortran_vec ();
  236.  
  237.       F77_XFCN (dgemv, DGEMV, ("N", nr, nc, 1.0, m.data (), ld,
  238.                    a.data (), 1, 0.0, y, 1, 1L));
  239.  
  240.       if (f77_exception_encountered)
  241.         (*current_liboctave_error_handler)
  242.           ("unrecoverable error in dgemv");
  243.     }
  244.     }
  245.  
  246.   return retval;
  247. }
  248.  
  249. // diagonal matrix by column vector -> column vector operations
  250.  
  251. ColumnVector
  252. operator * (const DiagMatrix& m, const ColumnVector& a)
  253. {
  254.   ColumnVector retval;
  255.  
  256.   int nr = m.rows ();
  257.   int nc = m.cols ();
  258.  
  259.   int a_len = a.length ();
  260.  
  261.   if (nc != a_len)
  262.     gripe_nonconformant ("operator *", nr, nc, a_len, 1);
  263.   else
  264.     {
  265.       if (nr == 0 || nc == 0)
  266.     retval.resize (nr, 0.0);
  267.       else
  268.     {
  269.       retval.resize (nr);
  270.  
  271.       for (int i = 0; i < a_len; i++)
  272.         retval.elem (i) = a.elem (i) * m.elem (i, i);
  273.  
  274.       for (int i = a_len; i < nr; i++)
  275.         retval.elem (i) = 0.0;
  276.     }
  277.     }
  278.  
  279.   return retval;
  280. }
  281.  
  282. // other operations
  283.  
  284. ColumnVector
  285. ColumnVector::map (d_d_Mapper f) const
  286. {
  287.   ColumnVector b (*this);
  288.   return b.apply (f);
  289. }
  290.  
  291. ColumnVector&
  292. ColumnVector::apply (d_d_Mapper f)
  293. {
  294.   double *d = fortran_vec (); // Ensures only one reference to my privates!
  295.  
  296.   for (int i = 0; i < length (); i++)
  297.     d[i] = f (d[i]);
  298.  
  299.   return *this;
  300. }
  301.  
  302. double
  303. ColumnVector::min (void) const
  304. {
  305.   int len = length ();
  306.   if (len == 0)
  307.     return 0.0;
  308.  
  309.   double res = elem (0);
  310.  
  311.   for (int i = 1; i < len; i++)
  312.     if (elem (i) < res)
  313.       res = elem (i);
  314.  
  315.   return res;
  316. }
  317.  
  318. double
  319. ColumnVector::max (void) const
  320. {
  321.   int len = length ();
  322.   if (len == 0)
  323.     return 0.0;
  324.  
  325.   double res = elem (0);
  326.  
  327.   for (int i = 1; i < len; i++)
  328.     if (elem (i) > res)
  329.       res = elem (i);
  330.  
  331.   return res;
  332. }
  333.  
  334. ostream&
  335. operator << (ostream& os, const ColumnVector& a)
  336. {
  337. //  int field_width = os.precision () + 7;
  338.   for (int i = 0; i < a.length (); i++)
  339.     os << /* setw (field_width) << */ a.elem (i) << "\n";
  340.   return os;
  341. }
  342.  
  343. istream&
  344. operator >> (istream& is, ColumnVector& a)
  345. {
  346.   int len = a.length();
  347.  
  348.   if (len < 1)
  349.     is.clear (ios::badbit);
  350.   else
  351.     {
  352.       double tmp;
  353.       for (int i = 0; i < len; i++)
  354.         {
  355.           is >> tmp;
  356.           if (is)
  357.             a.elem (i) = tmp;
  358.           else
  359.             break;
  360.         }
  361.     }
  362.   return is;
  363. }
  364.  
  365. /*
  366. ;;; Local Variables: ***
  367. ;;; mode: C++ ***
  368. ;;; End: ***
  369. */
  370.